Shocks in Vertically Oscillated Granular Layers 
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We study shock formation in vertically oscillated granular layers, using both molecular dynamics 
simulations and numerical solutions of continuum equations to Navier-Stokes order. A flat layer of 
grains is thrown up from an oscillating plate during each oscillation cycle and collides with the plate 
later in the cycle. The collisions produce layer compaction near the plate and a high temperature 
shock front that rapidly propagates upward through the layer. The shock is highly time-dependent, 
propagating through the layer in only a quarter of the cycle. We compare numerical solutions of the 
continuum equations to molecular dynamics simulations that assume binary, instantaneous collisions 
between frictionless, inelastic hard spheres. The two simulations yield results for the shock position, 
shape, and speed that agree well. An investigation of the effect of inelasticity shows that the shock 
velocity increases continuously with decreasing inelasticity; the elastic limit is not singular. 

PACS numbers: 47.40.Nm, 52.35.Tc, 45.70.Mg, 62.50.+P 
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I. INTRODUCTION 

A. Background 

A successful continuum description of granular media 
could exploit the powerful methods of standard fluid dy- 
namics to describe a variety of granular flow phenomena, 
but the validity of a continuum description of granular 
materials has been questioned and remains an open prob- 
lem Q . Several investigators have proposed equations of 
motion for rapid granular flow using continuum fields - 
number density n, velocity u, and granular temperature 
T (|r is the average kinetic energy due to random parti- 



cle motion) In one approach, the interaction be- 

tween grains is modeled with binary, inelastic hard sphere 
collision operators in kinetic theory to derive granular 
continuum equations to Euler W\, Navier-Stokes 6], and 
Burnett JJ order. However, attempts to use such contin- 
uum equations to predict properties of granular media 
have mainly been restricted to idealized systems that are 
difficult to produce in experiment ^2, 8, 9, 10| and to 
steady state or asymptotic time limits 0,11,(9, 11]. The 
viability of these continuum equations for experimentally 
realizable, three-dimensional (3D), time-dependent gran- 
ular systems has scarcely been examined. In this pa- 
per, we use 3D inelastic hard sphere molecular dynamics 
(MD) simulations as well as 3D simulations of contin- 
uum equations to Navier-Stokes order to investigate the 
strongly time-dependent properties of normal shocks in 
a vertically oscillated granular layer. 



Shock waves 



fields develop a mathematical discontinuity between the 
compressed region and the undisturbed region. In a vis- 
cous, conducting fluid, the fields vary continuously, but 
the front steepens until the shock thickness (the width 
of the region in which the fields change from the undis- 
turbed values to the compressed values) is on the order 
of a mean free path. 

When a supersonic fluid impinges perpendicularly onto 
a fiat, impenetrable plate, a shock forms normal to the 
plate. The normal component of the mean velocity (with 
respect to the plate) is zero at the plate, but the velocity 
of the undisturbed fiuid far from the plate is still super- 
sonic, fiowing towards the plate. The Mach number de- 
creases abruptly from greater than unity to near zero as 
the plate is approached. The density, temperature, and 
pressure also change rapidly in the region of the shock 
front. This normal shock propagates back upstream in a 
direction opposite that of the upstream flow velocity. 

This paper concerns normal shocks in vertically oscil- 
lating granular layers. A distinguishing feature of gran- 
ular materials is that granular flows reach su per sonic 
speeds under common laboratory conditions 13. IllL IT^. 
Hence shock formation, which is achieved only under ex- 
treme conditions in ordinary gases, is commonplace in 
granular materials. Shock propagation has been exten- 
sively studied in rarefied gases |il3t. tl4] , but the similar- 
ities and differences between shocks in ordinary fiuids 
and in gran ular media is the subject of ongoing research 
[ill llSjl . A recent laboratory and molecular dynam- 
ics study examined time-independent behavior of oblique 
shocks in granular flow between two closely spaced plates, 
and the results were compared with 2D simulations of 
Navier-Stokes order continuum equations [T]| . 



If the Mach number Ma (the ratio of the local mean 
fluid speed to the local speed of sound) is greater than 
unity at the point where a fluid encounters an obstacle, 
a compression wavefront is formed near the object and 
steepens to form a shock. In an ideal fluid with no viscos- 
ity or heat conduction, the wavefront steepens until the 



C. Model system 

The system consists of a layer of grains on an impen- 
etrable plate that oscillates sinusoidally in the direction 
of gravity with amplitude A and frequency /. The di- 
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mensionless peak acceleration is F = Att^ p A/ g, where g 
is the acceleration due to gravity. This simple oscillating 
system has been found to yield standing wave pattern for- 
mation T^, convection clustering Ts], steady-state 
flow fields far from the plate ^'^d shocks We 
examine shock formation and propagation in layers that 
are approximately 9 particles deep as poured except in 
Sec. Ill E, where layer depth is varied. Throughout this 
paper, the phrase "the layer" refers to the dense region in 
which the volume fraction of particles is greater than 4% 
of the random close-packed volume fraction, Vynax = 0.65. 

When the maximum acceleration of the plate is greater 
than g, the layer leaves the plate during each cycle. When 
the layer is off the plate, it is cooled by inelastic collisions 
between particles, while the particles are simultaneously 
accelerated towards the plate by gravity; this leads to a 
large mean velocity compared to the sound speed (See 
Sec. Ill C). We show in Sec. Ill that when the layer 
contacts the plate later in the cycle, the Mach number is 
much greater than unity with respect to the plate, i.e., 
the flow is supersonic. 

Shocks similar to those presented here were found for 
a wide range of F and /. Here we report results obtained 
from a particular set of parameter values: F = 3 and / = 
0.095\/.g/(j, where a is the diameter of a particle. This 
would correspond to a frequency of approximately 30 Hz 
for particles with a diameter of 0.1 mm. The particles 
have a coefficient of restitution e = 0.9, except in Sec. HI 
D, where e is varied to investigate the effect of inelasticity 
on shock propagation. The number of particles in the 
container per unit area of the bottom plate is lO/fi^, 
except in Sec. HI E, where the number of particles is 
varied to investigate shocks in layers of different depths. 
That is, the layer would be lOc deep if the particles were 
arranged in a simple cubic lattice. In actual packings 
seen experimentally, 10 particles/cr^ corresponds to an 
average depth of approximately ff = 9cr as poured [20l |. 
For consistency with previous investigations, we use H 
(the depth of the layer as poured) to express the depth 
of the layer throughout this paper. 

The values of F and / in this study were examined 
previously in experiments and MD studies of spatial pat- 
tern formation in shallow granular layers |20j. In the 
present simulations, pattern formation was intentionally 
suppressed by considering a container smaller than one 
full wavelength in either horizontal direction (with pe- 
riodic boundary conditions at the side walls); the ab- 
sence of patterns permits better visualization of the de- 
pendence of the fields on height z. 

We present 3D continuum and MD simulations of the 
oscillated granular system, and analyze the behavior of 
the fields n, T, and Ma. Section II describes the meth- 
ods used to simulate and analyze the oscillatory time- 
dependent flow. Section III reports our results, including 
an examination of varying inelasticity and layer depth, 
and Sec. IV presents our conclusions. 



II. METHODS 

A. Molecular dynamics simulation 

We use an inelastic hard sphere molecular dynamics 
simulation that was validated in previous studies |^ 
by comparison to experimentally observed standing wave 
patterns for varying control parameters. The collision 
model assumes instantaneous binary collisions. Linear 
momentum is conserved, but energy is dissipated through 
inelastic collisions. As in the continuum model (Sec. II 
B), the dissipation is assumed to be characterized by a 
single parameter, the normal coefficient of restitution e: 

e = -</^^n, (1) 

where u„ — (vi — V2) • (ri — r2)/|ri — r2|, and v's and 
r's are velocity and position vectors for a pair of col- 
liding particles before the collision (no superscript) and 
after the collision (with a superscript *). Particles are 
monodisperse, and |ri — and |r* — are identically 
a. In previous studies, the normal coefficient of resti- 
tution depended on w„ 20] . We investigate the effects 
of such a model in Sec. HI D. However, to be consistent 
with the continuum model, we use a velocity-independent 
coefficient of restitution for most of this study. We also 
neglect the surface friction between particles as well as 
between the particles and the plate. 

The particles are constrained between the bottom plate 
which oscillates sinusoidally between height z = and 
z — 2A, and a ceiling which is fixed a,t z — 200cr. In 
the two horizontal directions, the width is 20cr, and the 
boundary condition is periodic. Except in Sec. HI E, 
3937 particles were used; the number of particles was 
chosen to correspond to H = 9a and the total mass of the 
layer matches with that in the continuum simulations. 
In Sec. HI E, 1969 particles were used to correspond to 
H = 4.5(7, and 5906 particles were used to correspond 
to H — 13.5(7. The container has the same e as the 
particles, and the mass of the container is assumed to be 
infinitely large compared to that of the granular layer. 

Fields are calculated by dividing the container height 
into small bins of size a. The number density is the num- 
ber of particles found in a bin divided by the volume of 
the bin. The number density n from MD simulations 
is averaged over the same phase angle of the oscillatory 
state for ten cycles of the plate. The velocity v of each 
particle is used to calculate the mean velocity in each bin, 
u = (v), where the brackets represent the average over 
all the particles found in a bin at the same phase angle 
of the oscillatory state for ten cycles of the plate. Gran- 
ular temperature is then defined as T = i(|v — up). 
Throughout this paper, we use units such that the par- 
ticles have mass unity. 
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B. Continuum simulation 

We numerically solve the continuum equations of 
Navier-Stokes order proposed by Jenkins and Richman 
0. These equations were derived for a dense gas com- 
posed of frictionless (smooth), inelastic hard spheres 
by applying Grad's 13-moment method to the inelastic 
Enskog-Boltzmann equation. This model yields hydro- 
dynamic equations for number density (or equivalently, 
volume fraction v — ^na^), momentum, and granular 
temperature: 



which arises from the inelasticity of collisions between 
particles. The bulk viscosity is given by 



A 



n(jT^/'^G{v), 



the shear viscosity by 
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and the thermal conductivity by 
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where the components of the symmetrized velocity gra- 
dient tensor E are given by: Eij — i {djUi + diUj) . The 
components of the stress tensor P are given by the con- 
stitutive relation: 



P — 



(A - -fi)Ekk 



5,j + 2^jlE,, 



and the heat flux is calculated from Fourier's law: 
q = -kVT. 



(5) 



(6) 



To calculate the pressure, we use the equation of state 
and radial distribution function at contact proposed by 
Goldshtein et al. to include both dense gas and in- 
elastic effects: 



p = nT [l + 2{l + e)G{v)]. 



(7) 



where 



G{u) - ug^{v), (8) 
and the radial distribution function at contact, g^, is: 



5o(;^) 



(9) 



where Vmax = 0.65 is the 3D random close-packed volume 
fraction. 

Equations (|2I4|I differ from those for a compressible, 
dense gas of elastic particles by the energy loss term 



12 nT3/2 
ii(l„e^)^G(.), 
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In marked contrast to incompressible isothermal liquid 
layers, n and T for an oscillating granular layer vary 
by orders of magnitude in both space and time during 
each cycle of oscillation; hence the viscosities and ther- 
mal conductivity are very far from constant and must be 
calculated at each timestep. 

The equations are numerically integrated using a sec- 
ond order finite difference scheme on a uniform grid in 
3D and first order adaptive time stepping. For densities 
near close-packed, terms of order 1/G{v) are negligible 
in Eqs. (|12|) and (|13|) . For reasons of numerical stability, 
these terms were neglected and artificial dissipation '23| 
was introduced in the low-density region above and below 
the layer. In the bulk of the layer, where the shock forms 
and propagates, the density near the shock is high, and 
the artificial dissipation has little effect on the solution. 

As in the MD simulation, the granular fluid is con- 
tained between two impenetrable horizontal plates at the 
top and bottom of the container. We initially conducted 
the simulation in a cell with a ceiling that oscillates with 
the bottom plate at a height 200(t above the bottom 
plate. For H = 9cr, where most of these simulations were 
conducted, an average of less than one particle per oscil- 
lation reached a height of z > 80ct in the MD simulation. 
Changing the height of the cell from 200(7 to BOcr in the 
continuum simulation resulted in a signiflcant speedup in 
computation time, but did not change the shock profile. 
The fractional root mean square difference over one cy- 
cle between the shock location produced in the tall cell 
(200(T high) and the shorter cell (SOcr high) was less than 
1%. Therefore, for computational efficiency, the contin- 
uum simulations were conducted with a ceiling located 
at a height 80ct above the bottom plate, oscillating with 
the bottom plate. In MD simulations of shallow layers, 
H < 5, many particles reach heights of z > 80ct, so we 
use a cell height of 200(T in continuum simulations for 
H <5 (see Sec. Ill E). 

As in the MD simulations, we use periodic boundary 
conditions in the horizontal directions. Initial simula- 
tions were conducted with a container 20ct wide in the 
horizontal directions, as in the MD simulation. For rea- 
sons of computational efficiency, we reduced the horizon- 
tal width of the cell to llcr, which produced no significant 
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changes in the shock behavior. The width in the horizon- 
tal directions is llcr for continuum simulations through- 
out this paper. 

The continuum equations are solved in the reference 
frame of the container, with the container oscillations 
accounted for by a sinusoidal external forcing term in 
the momentum equation, Eq. |j2Jl. We use impenetrable 
boundary conditions at the plates: = in the frame 
of the plate. Boundary conditions for horizontal veloc- 
ities and temperature are set to match those found in 
the MD simulation. For most of the oscillation, their 
vertical derivatives at the plates were found to be negli- 
gible in the MD simulation, so for simplicity the bound- 
ary conditions used for the continuum simulation were: 
dux/dz = duy/dz = 0, and dT/dz = 0, although these 
derivatives are not identically zero in the MD simulation 
for the entire oscillation cycle of the plate. 



C. Shock tracking and shock width 

With each cycle of the plate, a normal shock wave is 
formed and propagates through the layer. This shock 
wave separates a compressed region near the plate from 
the undisturbed region which is still falling towards the 
plate. Both the undisturbed region and the compressed 
region vary in space and time. In Sec. Ill we find that 
when the shock is in the layer, granular temperature 
increases with height in the compressed region, drops 
quickly in the shock, and then increases again in the 
undisturbed region (Fig. The shock is identified in 
simulations as the region in which the temperature falls 
from its maximum at the top of the compressed region 
to the local minimum at the bottom of the undisturbed 
region. The shock width d is defined as the distance 
between the locations of this maximum and this mini- 
mum. For purposes of finding the propagation speed of 
the shock, the "shock location" is defined as the steepest 
region of temperature decrease, which is located near the 
center of the shock region (shaded regions in Fig. . 
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III. RESULTS 

A. Shock profile during formation and propagation 

Both simulations start with a slightly randomized flat 
layer near the plate which is then oscillated. A periodic 
state is reached after a transient state of several cycles 
of the plate. Although most of the particles move to- 
gether as a high-density layer, some particles are always 
found above and below the high-density region. We de- 
fine "the layer" as the region in which the volume fraction 
is greater than 4% of the random close-packed volume 
fraction Vmax- Below this volume fraction, there are too 
few particles for accurate averaging in the MD simula- 
tion. 



FIG. 1: Dimensionless temperature T/(gcj) and volume frac- 
tion V as functions of dimensionless height z/a at four times 
ft in the oscillation cycle. For each time, a picture from MD 
simulation is shown in the left column, with individual parti- 
cles color-coded according to temperature: high T in red, low 
T in blue, and the bottom plate of the container shaded solid 
gray. The right column shows horizontally averaged v (blue, 
labeled on the bottom graph) and T/{ga) (red, labeled on the 
top graph) as functions of z/a (ordinate) for the same four 
times. The plate is shown as a horizontal black solid line, re- 
sults from MD simulation are shown as dots, and continuum 
results are solid lines. At times when a shock is present in the 
graph, the shock region as located in continuum simulations 
(see Sec. II C) is shaded in gray. The width of the shock 
varies throughout the cycle, as discussed in Sec. Ill C. 
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During each oscillation, the shock forms and propa- 
gates at the same phase angle with respect to the plate's 
oscillation, so the dynamics of the system are fully char- 
acterized by the time interval between time ft — and 
one cycle later, ft = 1. We now investigate the formation 
and propagation of the shock by examining the behav- 
iors found in MD and continuum simulations for various 
times ft during the cycle (Fig. P). 

1. ft — 0: The layer falls towards the plate 

At ft = the container is at its minimum height. The 
particles, thrown off the plate in the previous cycle, are 
now falling towards the plate. The temperature is nearly 
zero for most of the material, increasing with height only 
in the low density region above the layer (Fig. . The 
effects of the artificial dissipation which is used in the 
continuum simulation are most pronounced at this time 
in the cycle, while the layer is dilated and above the plate. 
The shock formed in the previous cycle has propagated 
out through the top of the layer and is not present in this 
picture. The difference between the MD and continuum 
simulation at large heights has little effect on the shock 
behavior at later times in the cycle. 

2. ft — 0.15; Shock forms as the layer strikes plate 

The layer now begins to strike the rising plate. The 
layer compresses and particle velocities are randomized 
by collisions, increasing v and T drastically near the 
plate. This results in the formation of a shock region 
where z/ and T change rapidly from the compressed, high 
temperature region near the plate to the dilute, low tem- 
perature undisturbed region which is still falling towards 
the plate. At ft = 0.15 the shock is just beginning to 
form and is broad and not fully developed (see Sec. Ill 
C). During this time and throughout the formation and 
propagation of the shock, the shock profile from the con- 
tinuum model shows remarkable agreement with the pro- 
file found in the MD simulation. 



3. ft — 0.22; Shock propagates through the layer 

The layer continues to compress on the plate, result- 
ing in the propagation of the shock up through the layer. 
Here the shock is fully developed and has propagated 
through much of the layer, but there is still a significant 
number of particles falling towards the plate. The shock 
has steepened since ft = 0.15, with v and T changing 
from the upstream to the downstream values within a 
distance of about 2a. As will be discussed in Sec. Ill 
C, the width of the shock continues to change through- 
out the cycle as the shock propagates through parts of 
the layer with different densities. A discontinuity in the 
derivatives of v and T appears at the leading edge of the 



shock, showing the boundary of the undisturbed region. 
Collisions cause the layer to gradually cool behind the 
shock, creating a lower temperature near the plate. 



4- ft = 0.40; The layer begins to leave the plate 

The shock has propagated through the bulk of the layer 
and now enters the very dilute region above the layer. In 
this region the temperature becomes uncertain in the MD 
simulation because there are too few particles to use in 
averaging. Simultaneously, the plate is approaching its 
maximum height, and the layer begins to leave the plate 
as the downward acceleration of the plate becomes larger 
than g. The layer continues to cool, setting the stage for 
the next oscillation. 



B. Shock dynamics 

The formation and propagation of the shock is shown 
in Fig. 121 for three values of the restitution coefficient 
e. The effect of varying inelasticity will be discussed in 
Sec. Ill D. In each case, when the bottom of the falling 
layer hits the plate, the center of mass of the layer is still 
falling towards the plate, and the layer is compressed. At 
this time, a shock is formed near the plate. The shock 
propagates up through the center of mass and out into 
the low density region above the layer. Then the layer 
leaves the plate, dilates and cools during its free flight, 
and once again falls towards the plate. As the shock 
propagates through the layer, it goes through regions of 
differing v and T, and the shock velocity shows some 
slight variation with height, as seen in the curvature of 
the shock location (Fig.|21). 



C. Mach number and shock width 

A prerequisite for shock formation is that the local 
Mach number of the flow be greater than unity with re- 
spect to the object causing the disturbance. We calculate 
the speed of sound from MD and continuum simulations 
using a relation derived from the continuum equation of 
state Eq. [13: 

where x 1 + 2(1 + e)G{iy). 

As the layer falls towards the plate, it is accelerated by 
gravity and becomes supersonic in the reference frame of 
the plate. When the layer hits the plate, it is stopped by 
the impenetrable boundary, and the Mach number falls 
to nearly zero with respect to the plate (Fig. |3J. The 
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FIG. 2: Location of the shock (solid lines for continuum, triangles for MD) and the center of mass of the layer (dashed lines 
for continuum, circles for MD) as a function of time ft during one cycle of the plate (thick solid line) for particles with (a) 
e = 0.99, (b) e — 0.90, and (c) e — 0.85, starting from the same initial conditions at ft = —0.33. The plot is shaded according 
to the volume fraction from the continuum simulation, so that high volume fraction is dark and low volume fraction is light. 
The "top" and the "bottom" of the layer from MD (when the volume fraction drops to less than 4% of random close-packed) 
are shown as -|-'s and x's, respectively. The material below the shock is visibly compressed as compared to the region above 
the shock, as can be seen from the shading. 



undisturbed region is still falling towards the plate at su- 
personic speeds, and the shock is formed at the bound- 
ary between the compressed (low Ma) region and the 
undisturbed (high Ma) region. The Mach number in 
the undisturbed region is largest near the shock where 
the layer has been falling the furthest under gravity, and 
gradually decreases with increasing height. The Mach 
number is directly calculated from the values of v, u, 
and T, so the apparent increased discrepancy between 
the two simulations as compared to the discrepancy in 
previous figures is due to propagation of errors in these 
quantities in calculating Ma. 

While the shock is in the layer, the Mach number al- 
ways reaches its maximum at the lower boundary of the 
undisturbed region. Thus the maximum value of Ma 
at a given time represents the value at the leading edge 
of the shock. The flows in this system are hypersonic, 
with Ma > 10 just above the plate immediately be- 
fore the shock is produced (Fig. ^ . As the shock moves 
through the layer, the maximum Ma decreases since the 
fastest flow has been stopped by the shock, but the flow 
just above the shock remains supersonic until after the 
shock leaves the layer at time ft = 0.35, and Ma al- 
ways changes from supersonic to subsonic as the shock is 
crossed moving towards the plate. 

Since the density of the layer varies in space and time, 
the mean free path of the particles in the layer changes 
as a function of height and time. The mean free path ^ 
may be estimated as pi| : 



(15) 




n (jra^) 



4 6 
Ma, 10v 



FIG. 3: Mach number (Ma) (dashed line for continuum and 
triangles for MD) and rescaled volume fraction (lOi/) (solid 
line for continuum and circles for MD) as functions of z/a 
when ft — 0.22. At this time the shock is fully developed and 
is propagating through the layer (see Fig. . 



From ft — 0.12 to ft = 0.17 the shock is forming but is 
not fully developed. The dilute bottom of the layer hits 
the plate first at ft = 0.12, so that while the initial shock 
is forming, a higher density region is still falling toward 
the plate. As higher and higher density material hits the 
plate, the shock width d broadens (Fig. 0J as it forms 
from ft = 0.12 to ft = 0.17. After ft = 0.17, the shock 
steepens until the fully formed shock is about one mean 
free path in width, where it remains until the shock leaves 
the layer at ft = 0.35. After the shock leaves the layer, 
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FIG. 4: The maximum Mach number (Ma) (dashed line for 
continuum and triangles for MD) and the ratio of shock width 
d to mean free path ^ (solid line for contirmum and circles for 
MD) as functions of time ft. The shock begins forming at 
ft = 0.12 and leaves the layer at ft = 0.35. 



FIG. 5: Average dimensionless shock speed Vshocklyfg^ in 
the reference frame of the plate as a function of coefficient of 
restitution e. Vshock is calculated as the average speed of the 
shock from when the shock is formed until it leaves the layer. 



d remains approximately one mean free path, but the 
mean free path itself increases in the low density region 
above the layer, leading to a broadening of the shock 
front as the shock leaves the layer (compare ft ~ 0.40 to 
ft = 0.22 in Fig.UJ. 



D. Effect of inelasticity 

Previous studies have investigated shocks in elastic and 
inelastic materials in the asymptotic limit of infinite time 
and in inelastic materials at finite times ^3 for a 
piston moving into a uniform granular gas at a constant 
speed in the absence of gravity. We wish to investigate 
the effect of inelasticity on shock dynamics at finite times 
in the granular shaker system. However, the oscillatory 
state of the layer is not appropriate for comparing a gas of 
elastic hard spheres to one of inelastic hard spheres. An 
elastic gas, which has no mechanism for removing energy, 
will never reach an oscillatory state since the container 
injects energy into the layer at each cycle of the plate. In 
order to investigate the effect of inelasticity on shocks in 
granular materials, we slightly modify the system to look 
at the propagation of a single shock from identical initial 
conditions rather than looking at oscillatory behavior. 

We run simulations with different coefficients of resti- 
tution e, starting from identical initial conditions of a 
dilated layer falling towards the plate. These initial con- 
ditions are taken from the oscillatory state of the system 
for e = 0.9 using the same parameters as those used in 
previous sections. At time ft = —0.33, when the layer 
is near the peak of its flight above the plate, the coeffi- 
cient of restitution of the particles is suddenly increased 
or decreased. 

Changing e significantly changes properties of the 
shock and the layer following the next collision with the 



plate (Fig. [2J|. Since the initial conditions are the same 
for each case and the layer is in free fall, very few differ- 
ences result from the change in e until the layer collides 
with the plate. When the layer hits the plate it is com- 
pacted and many collisions occur between particles. The 
lower the coefficient of restitution, the more energy the 
particles lose through these collisions, and the more com- 
pact the layer becomes. In the case of high e, the layer 
rebounds and dilates quickly after the collision, while in 
the case of lower e, the layer remains very compact and 
leaves the plate almost as a solid body. The center of 
mass of the high e layer flies very high before beginning 
to fall back towards the plate, while the lower e value 
leads to a more shallow flight, with the center of mass 
beginning to fall earlier in the cycle. The shock trajec- 
tory exhibits curvature similar to that of the center of 
mass, remaining mostly straight for high e, but bending 
downward significantly for low e. 

After the collision with the plate, we track the propa- 
gation of the shock through the layer with a range of new 
values of e. Starting from the same initial conditions, the 
shock starts at the same time and height regardless of e. 
However, the velocity of the shock depends on e, as Fig.jS] 
shows. Although the shock speed changes somewhat as 
it propagates through the layer (see Fig.|21) , the average 
shock speed Vshock increases monotonically with e (see 
Fig. [SI . The special case of elastic particles appears to 
match with the limit e — + 1; thus there is no qualitative 
difference between elastic and inelastic gases. 

Throughout this paper, both MD and continuum sim- 
ulations were based on a collision model in which e was 
independent of impact velocity. A more realistic hard- 
sphere collision model could include a coefficient of resti- 
tution e = e{vn) that depends on the normal compo- 
nent of the relative colliding velocity of the particles, 
Vn, such that the collisions become elastic as Vn ap- 
proaches zero. To test the effect of such a velocity- 
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FIG. 6: Location of the shock (A for constant e = 0.9, x 
for velocity-dependent e = e{v„)) and the center of mass of 
the layer (o for constant e — 0.9, + for velocity-dependent 
restitution e = e («n)) as a function of time ft during one 
cycle of the plate (thick solid line) in the MD simulation. 



dependent restituition coefBcient, we ran MD simulations 
including e (ti„) of the same form as that used in studies 
of pattern formation in vertically oscillated layers [20j : 

e (vn) = 1 — 0.1 (vn/y/ga^)^ * for < y/ga, and e = 0.9 
otherwise. The velocity-dependence of the coefficient of 
restitution produces negligible changes in the shock pro- 
files and trajectories, and less than a 1% change in the 
average velocity of the shock (Fig. O . For the range of 
restitution values used in this paper (0.8 < e < 1.0), 
very few particles collide with low enough relative veloc- 
ity for the velocity-dependence to significantly affect the 
outcome. However, for velocity independent e < 0.8, the 
MD simulation was unsuccessful due to inelastic collapse, 
while velocity dependent e (u„) may be used to avoid in- 
elastic collapse with more inelastic particles. 



E. Shocks in layers of different depths 

Thus far, we have restricted our investigation of shocks 
in oscillated layers to layer depths oi H = 9a. We now 
investigate how changing H affects shock formation and 
propagation for particles with e = 0.9. 

For fluidized layers, increasing the layer depth in- 
creases the pressure compressing the layer on the plate 
during collision with the plate, resulting in tighter pack- 
ing of the particles near the plate. For 4(7 < iJ < 15a, 
increasing the layer depth increases the volume fraction 
of the layer near the plate, so that the trajectory of the 
center of mass of a layer of depth H = 13.5a is nearly 
identical to that for H = 9a, even though th e form er has 
1.5 times as many particles in the layer (Fig. 7(a) I. Both 
trajectories are actually lower than the center of mass 
trajectory for H = 4.5cr. 
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FIG. 7: (a) The height of the center of mass of the layer, 
and (b) the location of the shock as a function of time ft for 
various layer depths H. The trajectory of the center of mass 
(a) is shown for a complete cycle in the oscillatory state from 
the MD simulation with H = 4.5cr (□), H = 9a (•), and 
H = 13.5(7 (A). The center of mass location from continuum 
simulations agrees with MD to within a root mean square 
difference of 3% over one cycle in each of the three cases. 
The shock location (b) is shown for the fraction of the cycle 
in which the shock is in the layer, from the oscillatory state 
ol H = 4.5(7 (□ for MD, — • — for continuum simulation), 
H = 9a {• for MD, — for continuum), and H — 13.5a (A for 
MD, for continuum). 



For H = A.5a, the layer is almost entirely gaseous, so 
much so that many particles reach heights greater than 
z = 80(7 in each oscillation of the plate in MD simulations 
with a ceiling at z = 200(7. Therefore, we use a cell of 
height 200(7 for continuum simulations oi H = 4.5(7 (see 
Sec. II B). 

Changing layer depth also affects the formation and 
propagation of shocks in the layer. For more dilute layers, 
the shock forms earlier in the cycle, when the lower edge 
of the dilute layer collides with the plate. Thus, increas- 
ing layer depth causes the shock to form later in the cycle, 
at which point the shock then propagates faster through 
the higher densities found in deeper layers (Fig. 7(b) I. 
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For H = 13.5(7, the pressure on the layer during coUi- 
sion with the plate is so high that in some small regions 
near the plate, particles begin to form ordered lattices 
and the volume fraction increases beyond the random 
close-packed volume fraction. For deep layers H > 15<j, 
a significant portion of the layer is compressed into a 
solid lattice. The equation of state Eq. ([TJl does not al- 
low volume fractions v > Vmax, so continuum simulations 
cannot reproduce densities greater than random close- 
packed found in H > 15a. For H < 4(t, the particles 
do not move as a coherent layer, but rather act as a gas 
filling the container for the entire cycle of the plate. 

IV. CONCLUSIONS 

Shocks form in the vertically oscillated granular system 
with each collision of the layer with the plate. Proper- 
ties of these shocks captured in both MD and continuum 
simulations show good agreement even though the sys- 
tem exhibits a strong time dependence and large spatial 
gradients. Results from the continuum simulation agree 
with the MD results even for ranges where the validity of 
the continuum equations is in question, including when 
e is significantly less than one (e = 0.8) and when the 
volume fraction approaches close-packed. We also found 
that deeper layers exhibit denser packing near the plate 
and higher shock speeds than shallow layers. In addition. 



adding a simple velocity-dependence to the coefficient of 
restitution did not significantly change the behavior of 
the system. Finally, we found that increasing the coef- 
ficient of restitution causes no singular behavior in the 
shock velocity as e approaches unity. 

These results demonstrate that continuum equations 
can describe even a complicated three-dimensional, time- 
dependent granular system. These results also suggest 
directions for further research. A more realistic model 
of particle collisions could include friction between parti- 
cles. The simulations should be further tested by direct 
comparison to laboratory measurements on shocks. Fi- 
nally, an investigation of the role of a shock wave in trans- 
porting energy from the plate to the layer could provide 
useful information about the mechanism by which oscil- 
lations of the plate provide the energy to drive flows and 
patterns. 
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